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ABSTRACT 

We consider light propagation in an inhomogeneous irrotational dust universe with vanishing 
cosmological constant, with initial conditions as in standard linear perturbation theory. A non- 
perturbative approach to the dynamics of such a universe is combined with a distance formula 
based on the Sachs optical equations. Then a numerical study implies a redshift-distance 
relation that roughly agrees with observations. Interpreted in the standard homogeneous setup, 
our results would appear to imply the currently accepted values for the Hubble rate and the 
deceleration parameter; furthermore there is consistency with density perturbations at last 
scattering. The determination of these three quantities relies only on a single parameter related 
to a cutoff scale. Discrepancies with the existing literature are related to subtleties of higher 
order perturbation theory which make both the reliability of the present approach and the 
magnitude of perturbative effects beyond second order hard to assess. 
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1 Introduction 


The fact that cosmological observations do not conform to the predictions of Friedmann- 
Lemaitre-Robertson-Walker (FLRW) models with a vanishing cosmological constant A is usu¬ 
ally interpreted as an indication that A differs from zero. Clearly our actual universe deviates 
from the idealized FLRW cases by hosting inhomogeneities, and there have been many sug¬ 
gestions that the latter might have effects which would explain the data without requiring A; 
see e.g. Ref. [1] for an early proposal of this kind. The main challenge for any such claim is 
to explain why we perceive an accelerated expansion. Basically there are two possible routes 
as well as combinations of them. On the one hand the inhomogeneities might have an impact 
on the actual expansion of the universe (suitably defined in terms of the evolution of volumes 
of spatial regions). On the other hand there is the possibility that they affect light propaga¬ 
tion in a subtle way which modifies the usual distance-redshift relations. In the present work 
we are mainly concerned with the second scenario, which relies on the obvious yet important 
insight that almost every single piece of evidence on the evolution of the cosmos relies on the 
observation of photons with telescopes or other devices; Ref. [2] provides a particularly forceful 
presentation of this point. 

There is an extensive amount of literature on light propagation in the presence of inho¬ 
mogeneities; see e.g. Refs. [3, 4, 5, 6, 7, 8, 9, 10, 11] for a small subset. Typical ingredients 
include the use of the Sachs optical equations [12] from which a formula for the angular diame¬ 
ter distance dn can be derived, and approximations of the Dyer-Roeder type [13]. A somewhat 
different approach is pursued in Refs. [14, 15, 16, 17] and related papers, where a tailor-made 
coordinate system [18] is used. 

The present work will take the Sachs optical equations as a starting point, but will use 
them to analyse the evolution of the “structure distance” (cf. Weinberg [19]) ds — (1 + z)dA- 
The result, a second order ordinary differential equation, looks more complicated at first sight 
than the corresponding formula for d^, but it turns out that the two nontrivial coefficients have 
very simple interpretations: one of them is a local (and directed) expansion rate that agrees 
with the standard Hubble rate in the homogeneous case, and the other one is a quantity that 
vanishes in a spatially flat homogeneous geometry. These expressions (more precisely: their 
suitably defined expectation values) are then computed non-perturbatively in the framework 
of a recently introduced statistical model [20] whose only assumptions are an irrotational dust 
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approximation for the matter content and initial conditions consistent with linear perturbation 
theory with only Gaussian fluctuations. With the help of some approximations (but not of the 
Dyer-Roeder type) and the use of a computer program we find that in such a universe with A = 0 
there is a time t Q with the following properties. An observer at t Q will see redshift-distance pairs 
which, if interpreted with formulas that ignore the inhomogeneities, would indicate H{t 0 )t 0 ~ 1, 
a deceleration parameter q(t Q ) ~ —0.5, and density perturbations at a redshift of z ~ 1090 from 
t Q that agree with those assumed for dark matter at last scattering. In other words, such an 
observer sees what present day cosmologists see, despite living in a universe in which the 
cosmological constant vanishes. 

In the next section we derive a differential equation for the structure distance and discuss 
the meaning of its coefficients; furthermore we elucidate the relationship between local ex¬ 
pansion data along a lightlike geodesic and the inferences that a cosmologist who ignores the 
inhomogeneities would make. In Sec. 3 the coefficients are computed explicitly for the cases of 
homogeneous and irrotational dust universes. Sec. 4 contains a brief summary of the methods of 
Ref. [20] for a non-perturbative statistical treatment of an irrotational dust universe with initial 
conditions from linear perturbation theory. In Sec. 5 the “photon path average” is introduced: 
this is the concept that we use to estimate the overall effect of the changing environments that 
a photon experiences on the way from its source to an observer. Sec. 6 contains calculations 
up to second order in perturbation theory (we will see that they do not suffice to produce the 
relevant effects). In Sec. 7 we present the results of a numerical computation that transcends 
perturbation theory: we find quantities that are in rough agreement with today’s observations 
even though we assume A = 0. In the final section we briefly reiterate our hirelings and sum¬ 
marize the approximations that were made in deriving them. We also explain why some of the 
approximations are not as good as they originally appeared, thus leaving the question of the 
non-perturbative impact of inhomogeneities still open; this is the main modification compared 
to previous versions of the paper. 
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2 Sachs equations and distance formulas 


Let us start with a brief summary of the homogeneous case in order to provide some reference 
points for our subsequent generalization. A homogeneous universe is usually described with 
the help of a time-dependent scale factor a(t) in terms of which the Hnbblc expansion rate is 
defined as 


Hit) 


a(t) 
a{t ) ’ 


( 1 ) 


and the deceleration parameter as 


a a 



-1 — 1 - 1 . 
dt\H 


( 2 ) 


The redshift z of a photon emitted at time t and observed at time t a , with both the source and 
the observer at rest with respect to a comoving frame, is given by 


1 + z 


a(t 0 ) 
a (t) 


(3) 


which implies 

H(t)=~Hl + z). (4) 

In the case of vanishing spatial curvature several distance formulas can be summarized as 

i=(1 + 2)1 |^' (5) 

where we have to take A = —1 for the angular diameter distance cIa, and A = 1 for the luminosity 
distance d l- The resulting identity dgjdA = (1 + z) 2 actually holds in any pseudo-Riemannian 
geometry; this is known as Etherington’s theorem [21]. The simplest version of Eq. (5) occurs 
if we take d to be the geometric mean of ci^ and dg, 


ds — (1 + z)dn — (1 + z) 1 d,L, 


( 6 ) 


for which there exists a variety of names in the literature; we will follow Weinberg [19] who 
calls ds the “structure distance”. Then A = 0, and Eq. (5) implies 


H = 


dz 

ddg 


(7) 


and, with Eq. (4), 


ddg = —(1 + z)dt. 


( 8 ) 
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In the following we consider an arbitrary spacetime geometry. We want to analyse a light¬ 
like geodesic corresponding to the path of a photon emitted at x/ and observed at x%. With an 
affine parameter s and a corresponding tangent vector k 11 = dx^/ds the redshift z is determined 
in general by the formula 

= (U • fc) e 

(u-k) 0 ' 

where u e and u Q are the normalized tangent vectors to the worldlincs of the source and the 
observer, respectively. If we assume that we have a distinguished timelike coordinate t such 
that both the source and the observer have worldlines with normalized tangent vectors d/dt, 
and that s is normalized so that ds = dt at the observer, we get 


dt 

1 + z = — , 
ds 


d d 

Le ' ds =(1 + z) Jt 


( 10 ) 


(to be evaluated at the source, i.e. at t — t e \ the same holds for the following equations). We 
write ^ or use dots when we treat t as parametrizing the geodesic, and we denote the partial 
derivative by the spacetime coordinate t = x° as do or j/. 


The Sachs optical equations [12] (see [22] for a textbook derivation) are 

-^f + ^pt' + kptl 2 = 


da, 


opt 


ds 


+ 2 9 opt a opt = --R a p l > u ,e 0 ‘k p e ,t k v , 


(11) 

( 12 ) 


where 6 opt and cr op t are the expansion rate and the shear of the null bundle, respectively. In 
general the terms expansion rate and shear refer to the change in the size and the shape of a 
bundle of geodesics. Since we will later apply the same notions to worldlines of dust particles, 
we indicate with the subscript that we are referring to the optical quantities. Furthermore 
£ = £(i) + a/—I£( 2 ) where £(i), £( 2 ) are spacelike unit vectors orthogonal both to k and to the 
observer’s worldline; because of these properties the right-hand side of the second equation 
remains the same if the Riemann tensor R Q 8 I1V is replaced by the Weyl tensor C' Q/ 3 pi ,, and 
corresponding effects are often referred to as “Weyl focusing”. The angular diameter distance 
dA is determined by 

-^ln d A = 0 opt , (13) 

which can be used to reformulate the Sachs equations as 


d 2 d 1 

= -(K P t| 2 + - R a pk a k l3 )d A , 


^(o'opt d\) = ^R aj3lxv £ a k^£^k v d 
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(14) 

(15) 



We now want to transform Eq. (14) into an equation for the structure distance ds — (1 + z)dA 
as a function of time. By using Eq. (10) we find 

ds - [ln(l + z)\ d s + id s = 0 (16) 


with 

i = (1 + ^) -2 (|cr 0 pt| _> + ^ R a pk a k p ) - ^ M 1 + z )- ( 17 ) 

As we will demonstrate in Sec. 3, the quantity i actually vanishes for spatially flat homogeneous 
universes. In that case Eq. (16) is solved by 


dst - 


(1 + z)dt = 


— [ln(l + z) 


dz. 


(18) 


Even for % ^ 0 the introduction of ds$ is useful because we can simplify Eq. (16) by treating ds 
as a function of g^#, which results in 

d 2 ds - 


dd 'U 

with boundary conditions at ds$ = 0 given by 


ds = 0 , 


(1 + zf ds 


dxU 


dd 


= 1. 


'5# 


(19) 


( 20 ) 


There is no perfectly natural way of generalizing the concept of a Hubble rate to an inho¬ 
mogeneous universe. Two operational definitions of a “Hubble rate” associated with a specific 
point on a geodesic can be made as generalizations of Eq. (7): 

dz „ dz 


//,nf = 


H,= 


( 21 ) 


dds 11 dds# 

Both formulas reduce to the standard Hubble rate for the case of a homogeneous spatially flat 
universe. While H in f is essentially the quantity that is inferred from observations under the 
assumption of flat homogeneity, is the expansion at the source in the direction of the photon 
emission: by virtue of Eq. (18) we have 


H, = -±lnll + z), 


( 22 ) 


in perfect analogy with Eq. (4); also note that 77 jj is just the second coefficient in Eq. (16). 
With the help of Eqs. (19) and (20) we End 


tfinf 


dds 

dds% 


1 + 



—i 

(1 + ^) 2 


d s dds{ — 


1 - 



i 

(i-M 


ds dt'. 


(23) 
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This means that the two definitions of H coincide at the observer, Hfit 0 ) = H in f(t 0 ) = H a , 
and that for positive i observations tend to overestimate and for negative i to underestimate 
expansion rates in previous epochs; in particular, for sufficiently large negative i we can perceive 
acceleration even if it does not take place. 


As we have seen, someone who ignores the nonvanishing of i (in other words, any cosmol- 
ogist believing in the standard concordance model) would interpret H m f as “the Hubble rate”. 
Furthermore, from Eq. (8) such a person would (wrongly!) infer a time parameter f in f with 


df inf 


ckU 


(U 


-dt. 


1+z l+z 
In fact, H'mt and t - m f satisfy an analogue of Eqs. (4) and (22): 

dz , , dz d 

H'vai — — = ~{1+Z)- 


dd s dt []■](' 

Let us also introduce the deceleration parameters 


dt\ 


inf 


ln(l + z). 


d 


<?inf = 


dt 


inf 


//, 


inf 


dt [ II. 


1. 


(24) 


(25) 


(26) 


By using the chain rule, the definitions of the various quantities and Eq. (16) one can show 
that they are related via 


<7inf = % + i 


■ ^s(l + z) 


ds z 


(27) 


This demonstrates again that negative i can lead to the perception of acceleration even if it 
does not take place. 


We can summarize the results of this section in the following way. From the values of the 
pairs (ds, z) along a given lightlike geodesic, without taking into account the quantity i that 
encodes the effects of curvature and inhomogeneity, one would infer an expansion history along 
that geodesic in terms of quantities tint, H m f and g„ a f. The actual expansion history along that 
specific geodesic is encoded by t, H$ and q§. The two sets of quantities are related by Eqs. (23), 
(27) and 


Hi n { dds = Hi j dds$ = dz, 

(28) 

df i n f = Hftdt = — dln(l + z). 

(29) 


In reality we have at most a single data point (ds, z) for any observed direction, and we require 
a statistical analysis. As we will see, even and q$ (suitably averaged over photon paths) can 
become quite different from the corresponding results from volume averaging. 
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3 Homogeneous and irrotational dust universes 


While all of our results up to now are exact in an arbitrary pseudo-Riemannian geometry with 
a distinguished timelike coordinate, we assume in the following that the metric can be written, 
in the synchronous gauge, as 

ds 2 = g a pdx a dx l3 = —dt 2 + gij(t,x)dx l dx (30) 


this is true for any homogeneous spacetime as well as for irrotational dust, where the dust 
particles have constant space coordinates x l . We want to express our quantities in terms of the 
spatial 3-geometry with the time-dependent metric g VJ . To distinguish it from the spacetime 
geometry we adopt the convention that an expression with greek indices or at least one index 
of zero or a left superscript of (4 ) pertains to the 4-metric g a p, whereas any other quantity, in 
particular the Ricci scalar R = R], refers to g t] . The connection coefficients for the metric (30) 
vanish if two or three indices are 0, and the non-vanishing coefficients are 

bop' 7^)(}()rj , T i0j T ij0 — (){ | (jjj , ' dtp'A: k/ji'A:, (31) 

with the notation for d/dx° = d/dt and more generally d for d/dx 11 , so that 

^ = d 0 + x l di. (32) 


The expansion tensor 9) and the scalar expansion rate 6 are defined by 


1 , 


9 j = ^9 lk do g kj , 


0 = 0- = ^, 
Vs 


(33) 


and the shear is the traceless part of the expansion tensor, 


= d ) - 




(34) 


The Riemann tensor R a ^s can be expressed in terms of the expansion tensor and the Riemann 
tensor Rijki of the spatial metric g^: 


RoiOj 

R()ijk 

^Rijki 


-gikd o 0j — Oikd^j , 

@ij|k @ik\j 5 

Rijki @il@jk T QikQjl') 


(35) 

(36) 

(37) 


with 9ij = gikO k j and with the vertical strokes denoting covariant spatial derivatives. 
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We now want to specialize our analysis of photon paths to a metric of the type (30), with 
the assumption that both the source and the observer are comoving: x l e = const, x l 0 = const. 
Since T° = 1, the 0-component of the geodesic equation is 

d 2 t 1 . - , dx l dx 3 . . 

i? + 2 {d °^an7 = 0 (38) 

or, upon division by (1 + z) 2 and application of Eq. (10), 

it< 39 > 

As x ,J ' is light-like and x° = t, the spatial part x l must be a unit vector with respect to < 7 ^, 


gijX l x J = 1 , 


whereby the previous equation becomes 


— ln(l + z) = - <7 ijX 1 x^. 

dt y J 3 13 


Similarly we can transform the spatial component 


^ + + = o (42) 

ds z J ds ds J ds ds 

of the geodesic equation into 

x 1 + ^ x 1 — Gkix k x l x l + 2cr]id + Y l - k x 3 x k = 0. (43) 

o 

Upon using this, together with (32), in the derivative of Eq. (41), we find 
- ^ ln(l + z) = (d 0 + x l di) + (docrij + x k d k a ij )x l x : ’ - 2 a^-x 1 - a H x k x l x l + 2a l k x k + Y l kl x k x l )x 3 . 


Note that up to now we have never used the Einstein equations 


<4 R, — A^ g a g — 8nG]yT a p. 


Let us assume that the spatial part of the energy-momentum tensor is proportional to the 
metric, = gijT k / 3, and that T 0 i = 0. This holds not only in the homogeneous case but also 
in the general irrotational dust case, where T tJ = 0. Then Eq. (45) implies that the spacetime 
Ricci tensor R a g must be of the same type, ^Rij = g t] ^R k /3 and R 0i = 0, so that 

= R m (k«) 2 + ijy k‘k> = (1 + z) 2 (R„„ + 1 < 4 >if]5); 


(46) 



in the last step we have used k° = dx°/ds — 1 + z and g t] k l k 3 = k^ + (/c 0 ) 2 = (1 + z ) 2 . With 
the help of Eqs. (33) - (37) this results in 

^(1 + z)~ 2 R a pk a kP = - l -d Q e + | - a 2 . (47) 

The traceless spatial part of the Einstein equations amounts to 


d 0 a) + 9a) + r) = 0, 

which implies d^a^ = —9aij/3 + 2a k a k j — r l3 , where 

„ R 

r ij — ~ ~^9ij 


(48) 


(49) 


represents the traceless part of the spatial Ricci tensor. Using this after inserting Eqs. (44) and 
(47) into (17) we get 


i 


(1 + z) 2 1cropt | 2 + R/6 - a 2 + {-aijO - 2 a\a kj - r i3 + 2a ij a k ix k x l )x l x 3 
+x l d i 9 / 3 + x k {d k aij)x l x 3 — 2ai 3 Y l kl x k x l x 3 . (50) 


This result is still exact within the irrotational dust framework and also for any homogeneous 
cosmological model. In the latter case it reduces to i — R/6 = K/a 2 with K e {—1,0,1} so 
that i/( 1 + z) 2 = K/a 2 is constant; thereby Eqs. (18), (19) lead to the well known distance 
formulas that involve sin or sinh functions for K ^ 0. 

Let us also note that the equation (15) for the optical shear is determined by 

= (1 + z?(\eo„ -+ 2r« + *WV« - -i (51) 

for any metric of the type (30), as one can ascertain by using similar methods. This expression 
vanishes for any homogeneous model. 
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4 Mass-weighted average 


If we knew the spatial metric g l3 in the vicinity of a given lightlike geodesic in an irrotational 
dust universe, we could now compute the redshift and the structure distance along that geodesic 
simply by solving Eqs. (41) and (16) with input from Eq. (50) (assuming we are also solving for 
cr 0 pt along the way). In practice we do not know the precise form of the metric and need to rely 
on a statistical model; in addition we have to make simplifications to keep the computations 
manageable. As we aim for results beyond perturbation theory, we choose the approach of 
Ref. [20] for our underlying statistical model. The present section is devoted to a brief summary 
of the relevant ideas and results. The central concept in this approach is the mass-weighted 
average [23] 

(X) mw (t) = — [ X(x,t)p(x,t)y/g(x,t) d 3 x (52) 


(X) mw (t) = - / X(x,t)p(x,t)y/g(x,t) d 3 x (52) 

m v J D 

of a scalar quantity X, where T> is a large domain (e.g. all of the visible universe), p{x , t) is the 
local mass density and 

m v = / p{x, t)yfg{xgt) d 3 x (53) 

Jv 

is the mass content of V. For the case of an irrotational dust universe, energy conservation 
implies 

(- p(x,t)y/g{x,t) ) = 0 (54) 

and therefore (<9oX) mw = <9o(X) mw . This makes it possible to evade the technical difficulties 
that arise with the more common volume average, where averaging and taking time derivatives 
do not commute. Nevertheless volume averages are easily computed within this approach as 


/ v\ _ (Xp ) mw 

( A )vol “— 

\f J /mw 


(Xa 3 ), 


here a is the local scale factor defined as 


a(t, x ) = 


P 

P(t,x) 


where p is an arbitrary fixed mass. Then the dust expansion rate can be expressed as 


9(t,x ) = — 


d 0 p(t,x) d 0 a(t, x) 


P(t,x ) 


a(t, x) 


and a set of rescaled quantities 


p = a 3 p, a) = a 3 a), 


R = a 2 R , r) = a 2 r !■ 
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obeys the evolution equations 


where 



d op = 0, d 0 dj = -ard 0 R = -2a 3 <J l j r j i , 

(59) 

=a- 3 (~l*i d + la + l^a) + 

(60) 

a = la + a - lad - a- 

(61) 


The initial values for these evolution equations can be found by comparison with linear per¬ 
turbation theory: upon neglecting vector, tensor and decaying scalar modes the space metric 
g\ ] ; (t, x ) at early times can be expressed in terms of a single time-independent scalar Gaussian 

random function C(x) as 


9if PT \t,x) = 4dsW 


5ij + 


y 1 3 




+ t^didjC{x] 


(62) 


here a^ds = const x t 2/ ' ! is the standard EdS (Einstein-de Sitter, i.e. flat matter-only FLRW) 
scale factor. By comparing with section 5.3 of Ref. [19] one finds that this metric is equivalent 
to a Newtonian gauge metric with <h = T = —C/3. It turns out that the initial conditions for 
our evolution equations are 


lim 4 

*-*•0 t3 

= (6nG N p) 1/3 , 

(63) 

d in (x) 

= 0, 

(64) 


20 2 

(65) 

Rini.^') 

= - — (6nG N p)^S(x), 

Bn )'j(x) 

= -^-(67TG N p)h lk s kj (x), 
y 

(66) 


where S and s k j are the trace and traceless parts of the matrix 

didjCix) = Sijix ) = Sij(x) + (67) 


of second derivatives of the function C(x). In this setup it can be shown that 


R(t ) = R m + 2 a 4 (f) o(t) + - 9(t)a 4 (f) a{t) dt, 


and that the evolution equation of the local scale factor a(x, t) is 


1-2 4 rt 


(d 0 a) =-nC N pa - -R in + -Aa -- / 9{t)a (t) a (t) dt. 


6 


9 


( 68 ) 


(69) 
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As long as one neglects the last term (a 2 Y kl 3 \ k ) in Eq. (60), the evolution in a given region 
will depend only on the initial conditions within that region; furthermore, if one chooses a 
coordinate system in which the symmetric matrix Sjj(x) is diagonal then r tJ and will be 
diagonal in that system at any time t. In this way it suffices to work with the probability 
distribution for the three eigenvalues of S tJ . As shown in Ref. [20], the assumption that C(x ) 
is a Gaussian random held suffices to compute this distribution explicitly in terms of a single 
dimensionful parameter which is related to the value of an integral that requires an ultraviolet 
cutoff. Then one can switch to dimensionless units by taking a specific value for this parameter. 
With the computationally convenient choice that was adopted in Ref. [20] and that will also be 
used here, one finds 

<S 2 > mw = 5, (. s ijSkl 8 ik 8i l ) mw = 10/3. (70) 

If one also chooses p such that 6ttGnP = 1 in the corresponding units then the perturbative 
series for a starts as 

2 Six) 4 S 2 (x) + 2s i j(x)s k i(x)8' lk: 8 jl 2 

a(x, t)=G + 3 -kG- ^ J - 1 2 + ..., 71 

6 84 

where we have neglected cubic and higher orders in perturbation theory. 

In the following we will develop the theory further in terms of the dimensionless quantities 
that we have introduced here. This leads to unique results (up to ambiguities in approxima¬ 
tions), with the only free parameter coming from the reintroduction of a dimensionful scale 
once we start comparing our results with physical quantities. 


5 Photon path average 

Finally we want to connect the distance formula (16), which relies on the values of the quantities 
H$ = — [ln(l + z)\ and i along a photon path, with the model of Ref. [20] as summarized 
above. We propose to do the following. We replace the right-hand sides of Eqs. (41) and 
(50) by suitable expectation values which we will denote by ( ... ) pp , where the subscript 
stands for “photon path”. The idea is that (X) pp (t) should be the average of X over all 
spatial positions x occupied by a photon of a given type (e.g. siipernova or CMB) at the 
time t, as well as all directions v of propagation of such a photon. A complete realization 
of this concept would automatically guarantee consistency with correct ensemble and angular 
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averaging. While the approximation we will make at the beginning of the next paragraph 
leads to a mild angular deviation, statistical isotropy and homogeneity will be manifest in all 
our computations. Every photon path corresponds to a random walk in the probability space 
determined by the six entries of Sij (or, alternatively, three eigenvalues and three direction 
components). Then X = (X) pp + AX with (AX) pp = 0, and by the linearity of Eq. (16) 
the contribution of AX gets small if a photon probes different regions of the probability space 
within a short time. 


Every photon path corresponds to a curve C in x-space (the M 3 parametrized by the spatial 
coordinates x 1 , x 2 : x 3 ) that ends at x Q . In the flat homogeneus case these curves are just 
straight lines. If the shapes of these curves were not altered by the presence of inhomogeneities, 
then our model would tell us how the basic parameters are distributed with respect to the 
euclidean metric dl 2 = 8 l] dx l dx^ along such a curve. We will make the simple approximation 
of assuming the same distribution even in the general case. As a next step we want to move 
on to a description that is based on physical time rather than euclidean length. We denote by 


dx 1 . j dt 

~dT =X Jl 


(72) 


the tangent vector to C normalized to euclidean unit length, i.e. 8ijV l vi = 1. Upon taking the 
g-norrn gijV l vi of v and using Eq. (40) we find 


dt = sjgijV l vi dl , 


(73) 


which reflects the fact that the photon flight time is proportional to the traversed distance as 
measured with the physical metric g. For any path segment of length dl we average over the 
three basic parameters of the model (indicated by ( ... ) mw ) and over all directions v, and 
weight by the time dt = sjgijV l id dl spent in such a segment. This results in 


Wpp = 


(Ss* x V9ij‘U i v’cPv) mw 
(Is 2 V9ijV i v i d?v) mw 


(74) 


where the integrations are taken over the unit sphere S 2 = {v : 8 l jV l v 3 = 1} in tangent space; 
if A" depends on x l explicitly, we make use of 


which follows from Eqs. (72), (73). 


x = 


VUiji'd'-' 


(75) 
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Our aim is to compute (X) pp for the nontrivial coefficients in Eq. (16), i.e. for the cases 
X = — [ln(l + z)]' and X = i. To this end we require integrals over S 2 of expressions that 
are polynomials in the v l except for the occurrence of factors of \Jg %3 v l vi. Since exact results 
would involve elliptic functions we work in a basis in which the metric is diagonal and write 


with 


Then 


9 = 


9ij 9^ij T 7ij) 

911 + 922 + .9.33 


, 7ll + 722 + 733 — 0. 


(v / 9i,l' , v) = f )Jg(l +7i,W)^ = 9 A/2 (! + J + ...) 


(76) 

(77) 

(78) 


on the sphere 5ijV l v J = 1. For each term in this expansion we require only integrals of polyno¬ 
mials in the v l , such as 


f S 2 (v l ) 2n d 2 v = 47 t /(2 n + 1), f s2 (v 1 ) 2 (v 2 ) 2 d 2 v = 47t/15, 

/ S 2 (^ 1 ) 4 (^ 2 ) 2 d 2 u = 47 t/ 35 , f s2 (v 1 ) 2 (v 2 ) 2 (v 3 ) 2 d 2 v = 47r/105. 


(79) 

(80) 


From now on we simply omit any terms that are of quadratic or higher order in the 7 ^. While 
this may look excessively crude, one can check that even in the extremal cases of one or two 
vanishing eigenvalues the error is at most around 15%. For the integral in the denominator of 
(74) this gives, upon using (77), 


>s 2 


a / gijV l v^d 2 v ~ 47 i^fg. 


(81) 


According to Eq. (41), — [ln(l + z)] = 9/3 + (7 t f£ l x ] . Since 6 has no direction dependence, 


e 


9 


9 


is 2 


x\/ 9 ijV ivj d 2 v = - a/ g ij v i v : >d 2 v « 4vrv^- 


's 2 


(82) 


In evaluating the second term we use the fact that is diagonal in the same coordinate system 
in which g l3 is: 


rh. 


ijx'&y/gijvivi = = Vg ( + 7«)(u*) 2 ) ( 1 - ^ + •••)• (83) 


v i=l 


z=l 


Upon restricting this to terms linear in 7 ^ and using the formulas (79) and (77) we get 


is 2 


16 

(Jijtfx 3 ^Jg ij v i vid 2 v « — 7r v / I(°'i 1 7ii + ^722 + <73733)- 


(84) 
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Combining our results gives 


(—[ln(l + -)]') PP « 


(v^(56> + 4^711 + 4q-|7 22 + 4o-f7 33 )) mw 

15 ( yfo) mw 


(85) 


Next we turn our attention to (i) pp . Since no direction is singled out, the expressions in the 
second line of Eq. (50), which are all odd under x l —» — x% do not contribute after averaging. 
The optical shear cr opt is determined by Eq. (15). The behaviour for small t Q — t is easily found 
to be (7 0 pt ~ ^.{t 0 — t)R a pi ll ,£ a kP£L i k u , i.e. well-behaved and vanishing in the limit t —> t Q . Under 
a 90° rotation —> £( 2 ), £( 2 ) —> —£( 1 ) the right-hand side of Eq. (15) changes sign, hence 

its photon path average vanishes and the behaviour of a opt resembles a random walk around 
zero. Near t = 0 we can use the results of linear perturbation theory as presented in Sec. 4 to 
find that the right-hand side of Eq. (51) behaves like f -8 / 3 , hence that of Eq. (15) like f -4 / 3 . 
Naively this would result in <j opt ~ t~ l and a contribution of type f -2 / 3 to Eq. (50), which is 
the same power as the leading (second order) behaviour of the other terms, as we will shortly 
see; because of the random walk nature it will however be suppressed. In the following we will 
neglect the term (1 + z) ~ 2 1 cr op t | 2 hi Eq. (50), but keep in mind that i will receive a moderate 
positive correction for intermediate redshift values; in particular we should remember that this 
makes our results more reliable for smaller than for larger redshifts. According to Eq. (68), 

8 ft 

R = a~ 2 R m + 2cr 2 + ~a~ 2 / 9(i)a 2 cr 2 di, ( 86 ) 

3 Jo 

where R,\ n = lini/^ 0 « 2 R- The contribution of (—&ij6 — rij)x l xi can be treated like that of 
GijX l xi before, resulting in 


f ■ 16 

/ (— Gift - rij)x l x J y/gijvhjJdJv « + rj)7n + 

Js 2 ih 

With slightly more work we also find 

[ -2 G^GkjiJx 3 / gijv'yidFv « -^vrv^[(a}) 2 (5 + 4y n ) + ...] 
Js 2 -Lh 

and 
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/ 2 G ij Gkix l x :) x k x l yJg ij v i v^d 2 v PS —7T V /^[(o-}) 2 (7 + 8711) + ...]. 
Js 2 

Putting the pieces together we obtain 


(87) 

( 88 ) 

(89) 




4 

9a 2 



6 {t)a 2 G 2 di 



4 

105 


[(7 g \6 + 7 r\ + 6 (<t}) 2 ) 7 h + . 


( 90 ) 
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Our formulas rely explicitly on the spatial metric g tJ in the diagonal basis. To obtain it from 
the quantities whose evolution is studied in Sec. 4 we use 


11 9 

2 ^o hi <7n = -g u d 0 gn = 9\ = - + a\ = (In a)' + a\ 


gn(t) = const x a 2 x exp ( 2 / a\(t)dt ) , 


which implies 


with analogous expressions for g 22 and g 33. Comparison with Eq. (62) shows that the constant 
must be the same in each case, and that setting it to 1 corresponds to a normalization where 

(o 2 ) = UpLRW 


6 Perturbative results 


Before proceeding to the results of a non-perturbative numerical computation, let us first assume 
that we are still so close to the EdS case that in most regions perturbation theory provides a 
good approximation. We work with the dimensionless quantities described at the end of Sec. 4. 
Again our first goal is the photon path average of the right-hand side of Eq. (41). From Eq. (71) 
we find (to the same accuracy as there) 


_ 1 / S(x) 2 13S 2 (x) + 12s ij (x)s k i(x)S lk S :)l 4 


9(x,t) =2 1 IH- -^-t 


1 3 + ... 


The approximation (81) is valid at linear order, and with Eq. (92) and the fact that crj is 
traceless we get 

/ := / VgijV l vid 2 v = 47ra + 0(2); (94) 

Js 2 

0(n) means an expression of n th or higher order in perturbation theory. Since the perturbative 
expansions 6 = 6^ + 6^ + 6 ^ + 0(3) and I = + I (]) + J (2 ^ + 0(3) have deterministic 

leading terms (i.e., 6 •A = (9^) mw and J (0) = (J®) mw ) and first order terms whose expectation 
values vanish (i.e., (9 lyl ' , ) mw = 0 and (l^) m w = 0), we get 

- <W> " W = 9 (01 + («< 2 » + n(C) mw + 0(3); 


Wpp = 


note that J (2) has dropped out at quadratic order so that Eqs. (70), (71), (93) and (94) suffice 
for computing 

<0>pp ~ 2t_1 - ( 96 ) 
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to the same order as a and 9 before. The approximation (84) implies 

4 

{aijX l x J ) pp « + 0-2722 + o-|7 33 ) mw 

at leading (second) order. This can be evaluated via 


a 1 = —a 


3.-4 m 


„ 1 _i 

art at ~ —t s r ps -4 3 Sll 

Jo 1 5 1 3 

(here and in the following equation we only consider leading orders), 


7ll = — - 1 « e 2 Io a l d * - 1 ^ 2 
9 


a\dt ~ t 3 Sn 


and Eq. (70); the result is 


{(TijX i X J )p P « —4T 


(97) 

(98) 

(99) 

( 100 ) 


Combining this with Eq. (96) we obtain 

Wp P fti ^ rl + ^. ( 101 ) 

where the approximation again neglects terms of cubic or higher order in perturbation theory. 


In order to compute (i) pp up to second order in perturbation theory we require the mass- 
weighted average of Eq. (90). We begin with 


/ /— Fi'm \ / in i 

( v9 o - )mw ~ (-) 


o /mw ' - \ /mw 

or a 


-f H(S(1 _ i(! S )) mw = )!)(S 2 }„. = » (102) 

where the approximations neglect contributions of third or higher order in perturbation theory; 
the linear term has dropped out upon averaging. All other expressions in Eq. (90) are explicitly 
of quadratic or higher order: with Eq. (98) we find 

1 


cr^-t-H^+ ...), 

\ [ 9a 2 a 2 di « -A _ 5( s 2 +...), 

a Jo 6 


and Eq. (99) together with 


r! = a 2 f} 


5 _ 4 

_ 9 r ' Sl1 


implies 


( a iQ + r i)7n + • • • ~ (g — g )t 3 4a (s^ + ...) — -4 s (s^ + ...). 


Combining all contributions and using Eq. (70) we arrive at 


(<) 


pp 


1 50 .4 1 22 1 4 1 10\ _ 2 5 _ 2 

-x-(- (— x-x-x - — 4 3 = —t 3 . 

6 27 v 9 6 15 18 15 9 J 3 J 27 


(103) 

(104) 

(105) 

(106) 

(107) 
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7 Non-perturbative results 


In this section we present the results of numerical computations performed with GNU octave 
[25]. We used the Euler method with logarithmic time steps to solve the evolution equations 
(59) and (69). We assumed, however, constant r = f in instead of using Eq. (60), for the following 
reasons: the last term in that equation describes wavelike perturbations which probably play 
no role and cannot be described directly within the present model, and the other terms have 
extremely little impact on overall results (at least when volume evolutions are studied, see 
Fig. 12 of Ref. [20] and note that the tiny deviations only occur for t > 1). This was done 
for a large set of initial conditions, and the resulting values for a, a, r and R were used to 
evaluate the formulas of Sec. 5, with an appropriate probability measure for each set of initital 
conditions. More algorithmic details can be found in the appendix of Ref. [20]. 

In regions that collapse, the treatment in terms of irrotational dust breaks down and it 
is necessary to give a prescription on how to proceed with them. We followed the standard 
assumption, as suggested by the virial theorem, that collapsing regions shrink to half of their 
maximal sizes; somewhat unrealistically we pretended that such regions contract according 
to the irrotational dust evolution equations until that size is reached. The collapsed regions 
themselves were then treated in two distinct ways: firstly, by keeping them and letting all 
quantities retain the values that they had in the last moment of collapse, and secondly by just 
removing them from the statistics. The second approach makes more sense since it is doubtful 
whether many of the observed photons would have passed through a collapsed region, and 
also because the strong anisotropies that can occur during collapse should not persist in the 
virialized regions; nevertheless it is useful to have the other approach as well in order to get 
an idea of how strongly our results depend on details of modelling. In order to check that our 
results do not come solely from collapsing regions, we also performed computations in which 
we excluded any region from the statistics as soon as it started to contract. We will refer to 
these approaches as scenarios 1/2/3, respectively. 

The starting point is a computation of the basic results of the averaging process. The time 
evolution of (y/g ) mw = {\/{gn + #22 + 1733 )/3) mw as computed according to Eq. (92), does not 
differ substantially from that of its EdS equivalent t 2 / 3 (the discrepancy is less than 15% for 
the scenarios and time intervals that we consider here). We present our further results mainly 
in the form of figures created by GNU octave [25]. In these figures we use a colour coding of 
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blue/cyan/green for scenarios 1/2/3, respectively, with dashed lines for the quantities H$, q$ and 
ds$ and solid lines for the other quantities corresponding to these scenarios; furthermore EdS 
values are indicated by red dash-dotted, volume average results by solid yellow, perturbative 
results by dotted magenta and ACDM reference values by black dotted lines. 

1.6 


1.4 


1.2 

X 

1 


0.8 


0.6 

0 0.5 1 1.5 2 

time[units of Sec. 4] 

Figure 1: Time evolution of Ht 



Fig. 1 displays Ht over the time t for various versions of the Hubble rate H. The dashed 
lines in blue (highest), cyan (second) and green (third) correspond to the results of the non- 
perturbative computations; more precisely, they give (H§) pv t as computed numerically via 
Eq. (85) for the scenarios 1, 2 and 3, respectively. The fourth line (dotted, magenta) corresponds 
to the perturbative result (101), the fifth (solid yellow) line to Ht as computed via volume av¬ 
eraging, and the final red dash-dotted line shows the constant EdS value of H Ed $t = 2/3. 


19 






The strong deviations from the homogeneous case are a consequence mainly of local anisotropy, 
by the following mechanism. Consider a region TZ characterized by some specific values of 0 and 
c Tij and pick a frame {e!,e 2 ,e 3 } in which ay,- is diagonal. Assume, without loss of generality, 
that an > <r 22 and that originally TZ had the same diameters along the corresponding directions 
ei, e 2 . Even though the overall volume expansion of TZ is determined by 6, it will expand faster 
along ei and more slowly along e 2 , so that after a while 7 Z will have a larger extension in the 
ei-direction than in the e 2 -direction. A photon traversing TZ along ei will not only experience 
a stronger redshift per unit of time spent in TZ than one moving along e 2 , but it will also spend 
more time in TZ. The corresponding weighting that favors directions with stronger expansion 
results in the effect that on average a photon traversing TZ experiences a higher redshift than 
the volume expansion of TZ would suggest. 



Figure 2: Time evolution of i^fg 
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In Fig. 2 the time evolution of is displayed for our three non-perturbative scenarios; 
to be precise, (i)pp(y/g)mw, he. the mass-weighted average of the right-hand side of Eq. (90) is 
shown. The sharply dropping blue line corresponds to the first scenario, the curved cyan line to 
the second one, and the mildly dropping green line to the third one. These results are contrasted 
with the perturbative result iyfg = 5/27 and the EdS value of iyjg = 0 as represented by the 
two horizontal lines (in dotted magenta and dash-dotted red, respectively). Here the differences 
between the perturbative and non-perturbative results are not only enormous in magnitude but 
also change the direction of the effect. Once again the main contributions come from terms 
involving indicators of local anisotropy such as Oij and , as the form of the defining equation 
(50) suggests. 
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Figure 3: Structure distance over time 


Fig. 3 differs from the previous ones by relying not only on t e — t but also on t a , the present 


21 



age of the universe expressed in the dimensionless units of Sec. 4. Here and elsewhere our choice 
was simply to take t a as the time at which Hyt = 1 (remember that This is 

suggested by the fact that it seems to be a very good approximation in the case of the ACDM 
model and also close to lower bounds coming from ages of globular clusters; in a more general 
analysis one should probably also allow for values of H 0 t 0 somewhat above 1. For our first 
scenario we find ~ 0.7 in this way. The three lines ending at that value show various versions 
of the structure distance as functions of t = t e G [0,t o ]: the solid blue line shows dg itself, 
the dashed blue line below corresponds to dg$, and the dash-dotted red line that “starts late” 
corresponds to an EdS universe with the same value of H a , which would have had a shorter 
lifetime up to now. The other two triplets of lines correspond in an analogous way to the second 
scenario, where t Q ~ 1.35, and to the third one with t 0 ~ 2.3. 
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Figure 4: Structure distance over ln(l + z) 
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For producing Fig. 4, a plot of various versions of the structure distance over ln(l + z), the 
result of Eq. (85) (as shown in Fig. 1) was integrated to get ln(l + z) as a function of t, and 
combined with the values for the structure distance as displayed in Fig. 3. Each scenario is 
represented by a triplet of lines starting with the same slope which is lowest for the first and 
highest for the third scenario; the colour and linestyle coding are the same as before. This 
plot shows that < dst < dsi with differences of roughly the same size; i.e. the effect of 

a proper treatment of the second coefficient — [ln(l + z)]' in Eq. (16) is of the same order of 
magnitude as that of a proper treatment of the third coefficient, the quantity i. 



Figure 5: Deceleration over time 


Fig. 5 displays various versions of the deceleration parameter over the time t. The colour 
coding is the same as in the previous plots. The dashed lines give and the solid blue, cyan and 
green lines represent gw; in each case the lines end at our choice for t Q . The black dotted lines 
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correspond to deceleration in the standard ACDM scenario with = 0.72, with t 0 identified 
with the present time. Again the straight red dash-dotted line represents the EdS scenario, 
where q = 1/2, and the yellow line which shows only a slight downward slope displays the values 
that one gets via volume averaging. Once again we see that the photon path prescription leads 
to strongly different results, with effects of roughly the same order coming from the more precise 
treatments of the two non-trivial coefficients in Eq. (16). 

While all the results presented so far refer to times and distances in terms of the math¬ 
ematically convenient but observationally meaningless units of Sec. 4, the following plot uses 
standard units of years and parsecs. 
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Figure 6: Structure distance over time 


Fig. 6 is identical to Fig. 3 except for the normalization and the inclusion of a reference 
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ACDM curve (again as a black dotted line). This figure shows that, with the correct scaling, 
the predictions of the three different scenarios actually differ less than it appeared originally. 
Somewhat surprisingly, ds$ is closer to the ACDM values than ds here; in particular our results 
for ds overestimate the distances for early emission times. We can make this discrepancy 
quite precise by computing the distance to the last scattering surface from which the cosmic 
microwave background stems (see also Ref. [27]). This is not completely straightforward because 
the stepwidth of our programs is not fine enough for handling the time t\ s of last scattering 
that corresponds to z — 1090. We have circumvented this obstacle by using a combination of 
our programs and linear perturbation theory to find t\ s , noting that the solution of Eq. (16) 
near t — 0 takes the form ds(t) = ds( 0 ) +d^t 1 ^ 3 + 0(t 2 / 3 ), checking that the numerical results 
for small t are very well fitted by the first two terms, and using them to get ds(t\ 8 ). Upon 
doing this and converting the result to standard units, we found ds(t\ s ) ~ 20.7/20.9/19.8 Gpc 
for scenarios 1/2/3, respectively. These numbers overestimate ds by almost 50% compared to 
Planck results [28] of 13.9 Gpc (see their Table 2 and use ds = r*/#*[Mpc]), which is the largest 
discrepancy from standard values that we found in the present work. There are two possible 
explanations. On the one hand, we have omitted the term (1 + z)~ 2 \cr opt \ 2 (related to Weyl 
focusing) in Eq. (50); cf. the discussion after Eq. (85). Inclusion of this term would make the 
shape of the function ds(z) flatter and therefore more similar to the ACDM reference curve. 
On the other hand it is not clear whether the angular distance as inferred from the Planck 
results really should be exactly the same one as that computed via the Sachs equations. The 
Planck results refer to finite physical distances at t = t\ s , whereas the Sachs equations refer to 
the intersection of the observer’s backward light cone with that timeslice. In the homogeneous 
case this intersection will be perfectly spherical, but in a realistic inhomogeneous universe 
it might be somewhat crumpled (more like the surface of an orange), and the distance that 
corresponds to a total length along that surface (which is what the Sachs equations compute) 
will be somewhat larger. 
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Figure 7: Distance (normalized to EdS values) over z 


Fig. 7 displays, like Fig. 4, distance over redshift, the changes being the normalization of 
the distance to EdS values, the narrower range of z- values, the use of z instead of ln(l + z), 
and the inclusion of the ACDM scenario and supernova data. Again the red dash-dotted line 
corresponds to an EdS universe, the black dotted one to a ACDM universe with Da = 0.72, the 
solid lines to the observed structure distances d s for our three scenarios, and the dashed lines 
to the values of ds$- The black crosses mark the 551 supernovae from the Union2.1 compilation 
[29] that have z < 1, as taken from the Supernova Cosmology project website [30]. Both the 
second and the third scenario perform much better than the EdS case; actually the ACDM 
curve lies between the second and third scenario for most of the redshift values shown in the 
plot, and the second one somewhat overestimates the deviation from EdS. The first scenario, in 
which collapsed regions are included with the values for [ln(l + z)\ and i that they had in the 
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last moment of collapse, overestimates these deviations even more strongly. This suggests that 
our model would be improved by introducing a smooth slowing of the collapse (as it happens in 
reality), with a corresponding smooth transition of [ln(l + z)]’ and i to zero. The fact that even 
our third scenario, in which we have suppressed the effects from contracting regions, deviates 
strongly (and in the right direction) from the EdS case demonstrates that such an improvement 
could not obliterate the total effect of our treatment of inhomogeneities. 

What have we seen up to now? Considering a universe with A = 0 and with distributions of 
geometric quantities that follow directly from initial conditions based on a Gaussian distribu¬ 
tion, and with photons that obey the Sachs optical equations, we have shown that the following 
facts hold: there is a time t Q such that an observer at that time sees a redshift-distance rela¬ 
tion remarkably similar to that predicted by the standard ACDM scenario, and if the observer 
analyses the data without taking into account the inhomogeneities, he will infer a Hubble rate 
H m { such that H ini t 0 = 1 and a deceleration parameter q ini « —0.5. 


We have already considered the time t\ s of last scattering in our discussion of Fig. 6. We can 
make a further, less ambiguous, statement on that era in the following manner. In our most 
realistic scenario (the second one), t\ s ~ 5.3 x 10 -5 in the dimensionless units of Sec. 4 (with 
t — 0 the instant at which the singularity would have occurred in a purely matter dominated 
universe). At this time linear perturbation theory is still perfectly valid so that we can compute 
density perturbations at last scattering with the help of formulas (71) and (70): 



Is 




-tl AS 


^ x (5.3 x 10“ 5 )i x 


V^5 « 1.6 x 10~ 3 . 


(108) 


These are the density perturbations for the total matter, which are dominated by the ones 
for dark matter. According to Eq. (2.6.30) of Ref. [19], the density perturbations of baryonic 
matter satisfy Ap B /p B = 3AT/T, where T is temperature; using the commonly cited value 
of 1CT 5 for the relative temperature fluctuations in the CMB we find that the total density 
perturbations are roughly 50 times as large as those for the baryonic matter. This fits very 
well with the fact that dark matter decouples from photons (hence clumps gravitionally) earlier 
than baryons. Similar values for the ratios of the baryonic versus total density perturbations 
are required for structure formation; see e.g. Fig. 1 of Ref. [26]. We can turn this argument 
around: from the density perturbations we see that the time of last scattering cannot have 
occurred significantly before the time t\ s ~ 5.3 x 10~ 5 that corresponds to t Q ~ 1.35. But 
then it is clear that the inhomogeneities will have a significant impact on inferred Hubble and 
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deceleration rates, so that the assumption that a homogeneous universe (with or without a 
cosmological constant) give correct predictions necessarily breaks down. Conversely, since we 
do not require a non-zero A to account for present observations the simplest assumption is to 
take A = 0. 

8 Discussion and outlook 

Let us start our discussion with a brief reiteration of our assumptions and conclusions. Con¬ 
sidering a universe that 

• is matter dominated and obeys the Einstein equations, 

• in its early stages was very close to being spatially flat and homogeneous, with only 
Gaussian perturbations, and 

• has vanishing cosmological constant, A = 0, 

we found that there is a time t Q such that observervations made at that time and interpreted 
with formulas appropriate to the homogeneous case, would suggest 

• an inferred Hubble rate H [nf such that H in{ t 0 ~ 1, 

• an inferred deceleration parameter of q- m f ~ —0.5, and 

• density perturbations at a redshift of 1090 that fit well with values required at last scat¬ 
tering to lead to structure formation. 

In other words, an observer at time t Q in such a universe sees essentially what present day 
cosmologists see, even though A vanishes. This is the consequence of a model that has only 
one parameter (the overall scale) which can be adjusted. Once this parameter has been fixed 
by any of the three quantities that were just mentioned (and thus t Q identified with the present 
age of the universe), the prediction for either of the other two provides a highly nontrivial test. 
Our methods have performed very well on both of them. 



In order to arrive at these results it is essential to consider the effects of inhomogeneities on 
light propagation (not just on the evolution of volumes), and to use a formalism that transcends 
perturbation theory. The main steps involve the derivation of the differential equation (16) for 
the structure distance ds = (1 + z)dA, and the computation of the two non-trivial coefficients 
= — [ln( 1 + z)] and i that occur in this equation. In the spatially flat homogeneous case 
is just the usual Hubble rate and i — 0; otherwise each of these coefficients contributes 
significantly, with effects of roughly the same magnitude, to the deviations in the values of ds, 
H- m f and q- m {. The main source of discrepancies from FLRW universes is the local anisotropy, as 
encoded in the dust shear a tJ and the traceless part r t] of the Ricci tensor, and not so much the 
inhomogeneity which manifests itself by variations of the expansion rate 6 and the spatial Ricci 
scalar R. While Eq. (16) is valid in an arbitrary geometry in which photons follow light-like 
geodesics, the subsequent computations required a number of approximations: 

• The matter was modeled as irrotational dust. While this is an excellent approximation 
during expansion, it would not permit stable structures such as galaxies and clusters as 
the results of collapse. Our way of treating this problem, by simply assuming that collapse 
holds at half the maximum size (or ignoring collapsing regions altogether), is certainly 
somewhat ambiguous. In particular, the differences between the three variants that we 
chose show that the results do depend on such details; at the same time our third scenario 
demonstrates that deviations from the homogeneous case do not stem exclusively from 
collapse. As we argued in the discussion of Fig. 7, a smoother transition to the virialized 
state in our framework would probably lead to even better agreement with observations. 

• We have replaced statistical quantities by their expectation values in order to arrive at 
a description in which distance can be seen as a function of redshift, as in homogeneous 
models (cf. the first paragraph of Sec. 5). From the set of supernova data it is clear that 
this is a gross oversimplification. 

• We assumed a distribution of photon paths in x-space (the space in which our matter is 
at rest, which starts out as being almost perfectly euclidean) that was the same as if the 
photons moved along straight lines in that space. 

• While exact evolution equations were used for the local scale factor a, the shear and 
the Ricci scalar R , the evolution of the traceless part r VJ of the Ricci tensor was simplified 
by ignoring the right-hand side of Eq. (60). 
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• In our analysis of expressions that arise upon taking photon path averages, we have 
neglected terms of quadratic or higher order in 7 ^ (a scaled version of the traceless part 
of the metric < 7 ^). 

• For reasons that we discussed after Eq. (85) we ignored Weyl focusing, i.e. the contribution 
of the optical shear <j opt . 

• For the numerical treatment the time axis and the probability distribution for the back¬ 
ground parameters were discretized. The resulting errors are, however, much smaller than 
those coming from the other approximations. 

Unfortunately the second and third item are not as harmless as they originally seemed. Upon 
replacing i and by their expectation values, we have introduced errors A i and AH$ which 
have vanishing expectation values and are of first order. These lead to errors Ad and A ln(l + z) 
of the same type. The transition from d(t) and ln(l + z)(t) to d(z) then generates products 
of errors which are of second order and nonvanishing expectation value. The approximation 
introduced in the second paragraph of Sec. 5 probably leads to similar problems, whereas our 
computations in Sec. 6 respect the essential terms at second order perturbation theory. 

A general n th order term is an n-fold product of C (or, equivalently, the Newtonian potential 
<F) or its derivatives, in such a way that typically the n th order term has a total of up to 2{n— 1 ) 
spatial derivatives more that the first order term (see Ref. [24] for a detailed discussion). While 
C itself is small, d 2 C can be large; for example, density perturbations are of this type. In 
particular, among the terms contributing to the redshift-distance relation at second order, 
the largest ones that we find are proportional to ((<9 2 C) 2 ). However, according to the two 
independent groups that have performed complete computations up to second order [ 8 , 10, 15, 
16], terms of that type cancel out completely and subleading terms give corrections of an order 
of magnitude of only around 10 ~ 4 . 

This can be explained in the following way. In our approach, using the synchronous gauge, 
the whole setup relies on expressing quantities in terms of the entries (or eigenvalues) of the 
matrix Sij = didjC ; to be precise, the n th order contribution to any of the quantities a, < 7 , R and 
f is homogeneous of degree n in S. Terms of this type also produce the dominant contribution 
to the deviation of the metric from the FLRW case. But terms of (schematically) type d 2n C n 
in gij give rise to terms of type d 2n+2 C n in the curvature, which must all cancel. This means 
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that the part of the spatial metric consisting of the highest derivatives is flat, implying that it 
is possible to reparameterize the spatial slices in such a way that the metric no longer contains 
the d 2n C n terms. Hence any approximation in the synchronous gauge that does not respect 
the precise structure of the d 2n C n terms introduces errors that are potentially larger than the 
physical effects from the inhomogeneities. Since our approach suffers from this problem, it 
does not provide a conclusive argument that the standard ACDM picture require modification. 
Nevertheless it is intriguing how well it appears to perform - after all, one would expect mere 
errors to result in random nonsense rather than something that closely resembles observations. 
Besides, standard perturbation theory cannot be trusted either: higher order terms are not 
smaller than first order terms [24], and the real universe features shell crossings and vorticity, 
which do not occur in a purely perturbative modelling of an irrotational dust universe, but 
whose effects are taken into account by the approach to virialization in the present framework. 
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